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Abstract 

We exactly rewrite the Z(2) lattice gauge theory with standard plaque- 
tte action as a random surface model equivalent to the untruncated set of 
its strong coupling graphs. We simulate such surfaces including Polyakov 
line defects that arc moved by worm type update steps. Our Monte Carlo 
algorithms for the graph ensemble are reasonably efficient but not free of 
critical slowing down. Polyakov line correlators can be measured in this 
approach with small relative errors that are independent of the separation. 
As a first application our results are confronted with effective string theory 
predictions. In addition, the excess free energy due to twisted boundary 
conditions becomes an easily accessible observable. Our numerical experi- 
ments are in three dimensions, but the method is expected to work in any 
dimension. 
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1 Introduction 



This paper continues our series where we try to build on and generahze an idea by 
Prokof 'ev and Svistunov [1] on alternate ways of formulating and simulating simple 
statistical systems and Euclidean lattice field theories. In our context, their idea 
has two aspects. One is that, instead of summing over the field configurations in 
the original form, one exactly rewrites the models as a finite or infinite sum over 
associated strong coupling or high (classical) temperature graphs. The second 
aspect for spin models is, that enlarging the class of such graphs from those of 
the partition function to a larger set associated with the two point correlation 
function makes simulations ('diagrammatic Monte Carlo') both simpler and much 
more efficient with respect to critical slowing down. In [2] it was shown that we may 
in addition exploit some freedom to design the enlarged ensemble to our advantage 
and thus achieve an excellent signal to noise ratio for interesting observables. It is 
this aspect that we will extend in this paper to Polyakov line correlators in Abelian 
gauge theories while critical slowing down is unfortunately not eliminated at the 
same time. 

Some remarks are in order here. The standard usage of strong coupling ex- 
pansions is to evaluate a truncated series for some suitable observables and then 
take the thermodynamic limit within this approximation. These expansions in 
powers of some /3 usually have finite radii of convergence related to singularities 
or phase transitions. In systems with a large but finite number of compact de- 
grees of freedom however, these expansions converge for many quantities like for 
instance the partition function itself. Therefore, for such systems - and we never 
simulate anything else - the untruncated set of strong coupling graphs furnishes 
an exact reformulation of the theory in which we can attempt to apply stochastic 
summation methods. For a given action or Hamiltonian the graphs derived from 
it reproduce finite lattice observables exactly and not just universal features. For 
spin models we thus encounter representations in terms of loops drawn on the 
lattice and the obvious generalization to gauge models studied here consists of 
surfaces. The admissible graphs are restricted by constraining rules, for instance 
on the number of lines that may touch at a site, and we sum over such constrained 
objects. Solving the constraints in terms of independent new variables would com- 
plete our reformulation to a duality transformation, which is however not done 
in [1]. What the reformulations have in common is that rather different observ- 
ables may be advantageously computed in one or the other. In models similar to 
Ising models duality [3j relates the strong coupling expansions in one model to the 
weak coupling expansion in another. It is trivial, but amusing to remark, that a 
diagrammatic Monte Carlo of the weak coupling (low temperature) expansion of 
an Ising model is nothing but a standard simulation, as spin-flip related pairs of 
configurations correspond to one naive weak coupling graph. In the systematic 
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weak coupling approximation one just truncates these contributions guided by the 
size of their Boltzmann weights. 

Our previous work in this series extended the method from the Ising model 
[2] to the O(A^), CP(A^ — 1) models and to fermions, see [3] for an overview and 
more references. For fermions an unsolved sign problem hampers simulations in 
more than two dimensions while in the other cases the technique is expected to be 
efficient in arbitrary dimension. Studies of 0{N) as a special case of certain loop 
graph models on particular lattices are also found in [5], [Ej. 

There have been other recent attempts to apply the 'worm' method to Abelian 
gauge models [7j, [8], [9], [10], [11]. The important new feature in the present study 
is in our mind however, that we succeed in generalizing the low noise estimators 
for fundamental correlations from spin to gauge systems. Also the role of twisted 
boundary conditions on the torus will be discussed which leads to interesting topo- 
logical universal finite size observables. 

As we get access to precise results for Polyakov line correlators, a comparison 
with low energy effective string models suggests itself as a first interesting applica- 
tion of our simulation technique. Such studies form an active research area and we 
here offer a demonstration of the usefulness of our method in the context. A very 
profound such study is beyond the scope of this publication. The field has a long 
history with a boost, as far as we see, after the papers by Liischer and Weisz |12j, 
[I3] which interestingly were also triggered by algorithmic improvement. The low 
energy description of large Wilson and Polyakov loops in confining gauge theories 
was pushed to the two loop level in the spirit of a general effective field theory 
approach reminiscent of chiral perturbation theory vis a vis QCD. It was found 
that to this order and in three dimensions and taking into account all symme- 
tries, no free low energy coupling constants enter. This investigation was later 
even extended to the three loop level in [14] so that we have remarkable absolute 
predictions to check with our low noise long distance results for the Polyakov loop 
correlator. 

For the Z(2) Ising gauge theory in three dimensions that we adopt here as a test 
case, there actually exists an enormous body of precise results in the literature, 
see [T3], [TB] for recent papers with further references. Here the duality of the 
gauge theory with the spin model is exploited and the simulations are conducted 
there. The method is described in detail in pT|. A Wilson loop in the gauge 
theory translates into a ratio of partition functions in the spin model differing by 
flipped bonds. Such a ratio is factorized into many ratios differing in single bonds, 
which are evaluated in separate simulations. While this is an intrinsically three 
dimensional method, the same is not true for our more direct approach, although 
we here test it in D = 3 as well. 

The paper is organized as follows. In Section 2. we introduce the model, its 
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boundary conditions, rewrite it as a surface ensemble and collect some formulae of 
the effective string description. This is followed by the development of our Monte 
Carlo method for the surface ensemble in Section 3. and a report of our numerical 
experiments in Section 4. We end on conclusions in Section 5. and an Appendix 
listing numerical data. 

2 Z(2) gauge theory 

In this section we setup our model on a hypercubic lattice in D Euclidean dimen- 
sions. We include the definition of fluctuating twisted boundary conditions. Then 
we introduce the one-to-one reformulation as a surface model, which generalizes 
the loop (re)formulation of spin models. 

2.1 Twisted boundary conditions 

We consider a gauge field a{x,fi) = cr^X^) = ±1 G Z(2) defined on the links 
of a D-dimensional hypercubic periodic lattice of extent in directions fi = 
0,1, . . . , D — 1. The standard Wilson action is defined on plaquettes (x, /i < z/) by 

— S[a, P] = (3 P (x, /i, u) a{x, /i)cr(x + /t, z/)(t(x + z>, fi)cr{x, v). (1) 

X,IjL<V 

The plaquette dependent background field P(x, /x, z/) G Z(2) will be useful later 
and may be first imagined to be unity until further notice. Generalized periodic 
boundary condition^ demand that the gauge field is periodic up to gauge trans- 
formation^ |18] 

a{x + L^P, fi) = a{x, fi)T^ (x) r^, (x + fi) (2) 

specified by fixed 'external' transition functions t,^ (x). Because the shifts form an 
Abelian group it is necessary that shifts x ^ x + ^ x + + Lx come with 
the same gauge transformation as the double shift in the opposite order, which 
however still leaves the possibility 

Ti,{x)tx{x + L^u) = ju\rx{x)Ti,{x + LxX), 7i.a = 7ai. e Z (2) . (3) 

The case •yi.x = — 1 is possible because the action of r^, on gauge fields is indepen- 
dent of the global sign of the transition function. Then we have twisted boundary 

^ A fruitful point of view here is to consider all fields on the infinite lattice with a finite subset 
of independent variables as all others are 'locked' by periodicity. 

^ Gauge invariant densities will thus be periodic in the usual sense. This generalizes antiperi- 
odic boundary conditions in the Ising model which leave Z(2) symmetric composites periodic. 
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conditions in the z/A plane while we call planes with '-)yx = +1 untwisted or just 
periodic below. 

If we gauge transform 

a {x, fi) ^ a {x, fi) p{x) p{x + ft) (4) 

in the ordinary sense, this implies the following change of the transition functions 

(x) (x) p (x) p{x + LiyU) , (5) 

they transform like parallel transporters on a set of 'super-lattices' with spacings 
Lfj_. The 7^1/ represent the gauge invariant content of the transition functions. For 
a given set of twists 7^,^ we now define reference transition functions 

4''Hx) = llil.xt'^'^'^ (6) 

where we round upwards to an integer ('ceil') in the exponent. Then the product 
T,^ (x) rlf^ (x) has trivial twist and can be gauged to unity. Hence we may assume 
that the transition functions have the form t^"'^ and then the gauge field is periodic 
except (possibly) for 

a{x + L^O, p) = "f^^ai^x, p) if p < v and = (modL^). (7) 

By a further change of variables, we may arrive at fully periodic cr(a;, p) again if 
we absorb signs into 

P,(x,/x,z/) = (7^.)'^-»'-'» (8) 

with periodic 5 symbols (period L^^Ly respectively). Thus for each twisted plane 
we have a D — 2 dimensional set ('stack') of negative plaquettes that represent 
a background flux or disorder ('vortices') concentrated on a line, sheet, ... for 
D = 3, 4, . . .. This step is analogous to absorbing antiperiodic boundary conditions 
into a background gauge field, see [IS] for example. 

2.2 Gauge theory as a surface model 

Next we introduce the partition function with current insertions 

= 2-^' ^e-^l'^'^-^l JJa(x,/i)^'("'^). (9) 

In this formula we sum over strictly periodic a{x,p) = ±1 independently at all 
^ x^ < and p, and j{x,p) G {0, 1} is a periodic external field. The number 
of independent sites, links, plaquettes and 3-cubes (for later use) is 

N, = l[L^, Ni = DN,, Np = Ni{D-l)/2, = Np{D - 2) /3. (10) 
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Performing a local gauge change of variables ^ with periodic p (x) one derives 

z,[j] = z,[j]i[p{xr^'^^''^ (11) 

X 

with the divergence 

d;j,{x)^J2l^,ix)+j,ix-p)], (12) 

being the sum over the 2D links surrounding x. As p can be arbitrary this shows 
that only divergence free currents (in the Z(2) sense), for which 

d;j^{x) = (mod 2) (13) 

holds at all sites, yield nonzero Z^[j]. In addition we may flip the gauge field cr^ on 
any D — 1 dimensional hyperplane x^ = z orthogonal to the /i-direction ('layers'). 
A similar argument as before yields the requirement that the layer sums 

^j(a;,/i) = (mod2) ^ ^ < (14) 

must be even for all layers, which are hence pierced by an even number of current 
quanta. Equivalently we may say that the link field j {x,p) must have vanishing 
Z(2) winding number with respect to all torus directions. 
Using 

e^" = cosh(/3) ^ = tanh(/3) (15) 

n=0,l 

for each plaquette, we introduce a field n{x,p, v) G {0, 1}. Then we may average 
over the original a{x,p) which leave behind constraints only. We arrive at 

Z,[J] = 5^tS^..<^"(-''^''^)<|., [n] 6[d;n,,+j^], Z,[j] = {cosh pf^ Z,[j]. (16) 

n 

The sign is given by 

where for each n (x, /x, v) configuration we have introduced the 'wrapping' numbers 
Wf,y= ^ n{x,p,u) (mod 2), w^^ e {0,1} ■ (18) 
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They are topological quantitie^ that count (modulo two) how many times the 
surfaces made of plaquettes with n (x, /i, i/) = 1 wind around planes and represent 
two dimensional generalizations of the winding numbers of Ising loops. To define 
the divergence of the plaquette field, we extencQ 

n(x, /i, I/) = n(x, I/, yu), n(x, /i, yu) = 0. (19) 

The constraint (mod 2) combines the 2{D — 1) plaquettes surrounding each link 



We note that consistently this constraint can only be satisfied by j that obey (13) 



and (14). The proof looks slightly funny, for example 

dtju = d^dln^, = 2j2d;d:n,, dlj, = 0(mod2). 



An easy to interpret case is a configuration j = 0. Then the plaquettes n(x, /i, u) = 
1 form a surface which may branch and consist of disconnected components. The 
zero divergence condition demands that they are closed and have no boundaries, 
each link is surrounded by an even number of surface elements. Boundaries arise 



if j is nonzero and, due to (13) and (14) they may be written as a superposition 



of contractable closed current loops. It will be helpful for the reader to work out 
the closely analogous but geometrically simpler loop formulation of the Ising spin 
model in this language. 

In analogy to the steps taken for spin models we now form the current ensemble 

^ = J2R''[mj] = 5^t^--<'^"(^'^''^)i?-'[5x.] (21) 

j n 

where the non-negative weight R^^[j] will be specified later and Z without sub- 
script stands for trivial twist Z = Z^=i. Note that we here include all wrapping 
numbers in the sum over n with no extra signs, which corresponds to trivial twist 
in all planes. 

Expectation values in this ensemble are given by 

mn])) = |$^OMt^--"(^-'^'^)i?-^[9X.]. (22) 



■^We define it here also for configurations with j 7^ 0. A topological meaning, independent of 
our choice of r^'''-', is given however for vacuum configurations only. 

^Note that for mod 2 additive variables 'symmetric' and 'antisymmetric' coincides. 
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In addition we define 'vacuum' expectation values 



mn]))o 



{{0[n]S[d;n,^])) 



(23) 



on the subset of defect-free configurations. Such quantities do not depend on the 
choice of R. It is now obvious that ratios of partition functions of different twist, 
which lead to interesting observables, are given by topological observables 



Wc shall later see that there are Monte Carlo algorithms that are crgodic only 
in the sector of trivial wrapping in some or even all planes. It is clear now that 
the corresponding ensembles can be considered as arising by dynamically averag- 
ing over twisted and untwisted boundary conditions for these planes (fluctuating 
boundary conditions). 

By differentiating Z{{()[d*n^^])) with respect to /3 wc obtain the relation be- 
tween the average plaqucttc E of the original theory and the total surface area in 
vacuum configurations 



Below we shall find that close to the critical point in D = 3 we have the rather 
high values ((n))o ~ 1/3 and E w 0.95. 

2.3 Polyakov loop ensemble 

The source j {x, fi) allows to place a large variety of defect configurations like arbi- 
trary Wilson loops. We here define however a highly restricted framework involving 
only two Polyakov lines in the 0-direction that are located at {? = {ui, . . . ,ud-i) 
and an analogous v. The corresponding conserved current is 




(24) 




(25) 




(26) 



For coinciding u — v it vanishes and there is no defect. 
We consider the ensemble 



Z = J]p-^(«-t?)z[/">^)] 



■u,v 



■{u,v) 



(27) 



n,u,v 
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Here the weier ht has been speciahzed to < p{u — v) = p{v — u) while 



R ^ [j] vanishes if j is not of the form (26) 



The two point function (in the original gauge theory) of the Polyakov loop 
operator 

Lo-l 

7l{x)= Yl (To (x) (28) 
x'o=0 

is now given by a ratio of expectation values of counters 

G (x) = (vr (x) vr(O)) = p (x) (29) 

if we normalize p(0) = 1. 

2.4 Contact with transfer matrices, effective string theory 

The Polyakov line correlator is the ratio of partition functions with and without 
two static charges and thus given by 

G(f) = w„^e-^"(^)^°, (30) 

where e"^"*^"^-* labels the corresponding eigenvalues of the transfer matrix in the 
0-direction with charges separated by x and integer weights Wn account for degen- 
eracies [13]. 

Alternatively we may consider a transfer matrix in a spatial direction, k = 1 
for example. Then the Polyakov line operator excites a flux state whose energy 
strongly depends on its length Lq (now a 'transverse' direction) and we may con- 
sider the correlation 

Civ) = J2 Gi^) = Yl K\'e-^"' {y « L,), (31) 

■'■\^l=y 

where we have projected to zero momentum for the k > 1 directions and f„ is a 
nontrivial matrix element in this case. 

In effective string theories one can compute the energies £"„ as an asymptotic 
expansion in I/Lq. On the basis of the Nambu-Goto action for example one 
predicts for D = 3 the relation [20] [13] 



^- = s'il-l), s='A --=^ (32) 

Ss J vr 71 
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for the ground state energy which is expected to be relevant for z, s — )■ oo, and a is 
the (zero temperature) string tension. Formula (32) implies asymptotic expansions 



and, expanding s in z~^, 

It is interesting to note that in general the Nambu Goto picture is not expected 
to hold to all orders in the asymptotic expansion, but the terms exhibited above 
have been uniquely derived for general effective low energy actions just restricted 
by symmetries [13], [T4] - 

In addition it is shown in [13] that, up to rotational symmetry breaking cutoff 
effects, the original correlation is given in terms of the En by the expansion 

G ix) = Y: K\'2r ^ ir.(^_3) (E„r) (35) 

with r = and the modified Bessel function K^,. Note that here an infinite 
volume is assumed. 

The estimator for C (y) in our extended ensemble (27) follows from (29) and 
is given by the simple observable 

C^(y) = ^^%#|^ (36) 

that simply follows from the statistics of the defect-line separations. Note that this 
result is independent of the choice of p which we shall hence optimize with regard 
to the numerical simulation of (27). In practice one will symmetrize over spatial 
directions if they are symmetric with respect to extent and boundary conditions. 



3 Simulation algorithms 

We now introduce several algorithms to simulate the surface representation of Z(2) 
lattice gauge theory. 
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3.1 Defect conserving update 

A very simple update step consists of the following sequence CF ('cube flip') of 
operation^ 

• take a 3-cube c = x, fi < u < X, 

• propose to flip all 6 plaquettes associated with this cube which is the set 

z/, A) = 

z/), (x + X,fi,u), {x,fi,X) , (x + i>,/i, A) , (x,z/,A) ,{x + ft, z/, A)| 
i.e. n{p) 1 — n (p) for all p 

• accept this proposal with the Metropolis probability min(l,g) with 

Such steps may be iterated either with randomly chosen c or as a systematic 
sweep covering each c once in some order. This can be shown to lead to an ergodic 
algorithm at fixed j (vacuum graphs with j = for example) and within a fixed 
wrapping number sector. We expect (and confirm) however dynamical exponents 
close to two for such local updates that preserve the vacuum graph constraints. 
This is similar to just flipping links around plaquettes in the strong couphng loop 
representation of Ising spin models. 

A natural generalization of the worm algorithm [T] would be to allow to 'open 
up' vacuum graphs and allow excursions to the enlarged state space of (allowed) 
j ^ (possibly traversing 'useful' configurations regarding observables of interest) 
and eventually returning to a vacuum graph. We have a long record of not really 
successful experiments of this type, see P] for experiments with U(l) gauge theory. 
In Z(2) we have simulated correct ergodic algorithms where we included a chemical 
potential per defect link j {x, /x) = 1 to control their number. We could find values 
that have led to ensembles of randomly shaped (typically irregular) Wilson loops 
where also vacuum configurations (re-)appeared at a reasonable rate. None of these 
attempts has so-far led however either to fast dynamics or to easily interpretable 
observables. Therefore we come up with the hybrid below which combines cube 
flips with Polyakov line pair defects only where we have to tolerate however some 
critical slowing down. 

'^Such steps for U(l) have been introduced in [?!] . 
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3.2 Polyakov line moving update 



We now describe worm-type updates referring to the ensemble (27). If (say) u 
is moved to a neighboring site u' in the spatial direction i the corresponding line 
defect moves and the whole 'ladder' of Lq plaquettes in the plane 'above' this 
spatial link is flipped if the move is accepted. Such a move clearly preserves the 
constraint. To specify the procedure in detail it is helpful to define the auxiliary 
spatial link field 

A;(x,z) = ^ (x + 20,0,i) e [0,Lo) . (37) 

Now we execute the following sequence PS ('Polyakov shift') of (standard worm) 
steps: 

• \l u = V holds, we randomly re-locate both together to a new site on the 
D — 1 dimensional lattice. 

• Pick one of the 2[D — \) spatial neighbors u' of u, such that u' = u ± i. 

• Accept the proposal with the probability (pre-tabulated, of course) 

=min[l,t[^°-2^("'^V(M-t/)/p(M'-t/)] (u' = u+tj (38) 

or 

= min [l, t[^«-2K«-*^0]p (u - v) Ip {v! - v)\ (u' = u -tj . (39) 

If this happens, all plaquettes in the ladder are flipped and u is changed to 
u'. In the case of rejection the previous configuration is kept. 

3.3 Complete update sequence 

In this subsection we specialize to D = 3 for simplicity. One complete iteration of 
algorithm Al consists of N^/Lq repetitions of the following steps 

• Apply CF to ALq cubes around the temporal line at u. We use a helical 
order and first consider the 4 cubes around the link {u, 0) with uq = 0, then 
at Mo = 1 and so on. 

• This is followed by rips applications of PS. In all our simulations given below 
we took rips = 16. 
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The cost of one iteration Al is of the order volume hke a conventional sweep. We 
here try to mimic the successful worm algorithm in so far as we attach our update 
moves to the defect, which now is a (straight) line instead of a point. A difference 
in the case at hand is that defect moves alone - within our restricted class as 
discussed before - are by far not sufficient for ergodicity. We have mimicked this 
conglomerate also in the Ising model by proposing flips of links around only those 
plaquettes of which a defect forms a corner. This works but offers no advantage 
in this case. 

In the algorithm just described the defects at u and v are not treated on the 
same footing - we only shift u - although, due to the relocation step in PS they can 



both reach all positions. This is in contrast to (27) where they play a symmetric 
role. Analogous options appear already for the point defects in spin models, where 
we found the asymmetric algorithm easier to program and equally efficient as a 
manifestly symmetric variant. We make no such attempt here. 

As described before only temporal defect line pairs are included. Their migra- 
tions over the torus allow to change wrapping numbers in the 0/c planes, but not in 
the purely spatial planes (only 12 in D = 3). Therefore Al simulates an ensemble 
where one dynamically sums over 712 = ±1 with the other twists trivial. 

A second version A2 of the algorithm results, if we allow a pair of Polyakov 
lines in any of the 3 directions, but always only one pair at a time. We denote by 
j{^J■,u,v) ^YiQ current with lines in the /i-direction with u^v locating these lines. In 



(27) an additional summation over /i is included. 



^' = E E (40) 

and is the D — 1 dimensional sublattice of sites x with = 0. In addition the 



formulae in subsection 3.2 have to be modified in the obvious way. 

The practical difference between A2 and Al is simply that in PS, when u = vis 
encountered, a new value /i' G {0, 1, 2} is proposed together with a random position 
of u' = v' in V^'. The proposal is accepted with the probability min (1, L^//L^). 



This results in an ensemble as in (21) where all wrapping numbers contribute. 
If all extensions are the same the acceptance step is trivial and the resulting 
histograms of occurring u — v values can be combined. If this is not the case, 
they are collected separately and yield simultaneous results about Polyakov loop 
correlations in several geometries. 

3.4 Rejection free Polyakov shifts 

A potential problem with algorithm Al (and also A2) can be small acceptance 
rates for the shifts in PS, if Lq is large. After all we propose nonlocal albeit only 
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one-dimensional changes to the configuration. Simulating L'^ lattices at the critical 
point, we indeed find acceptances falling with L which reach values of about 1.6% 
only for L = 64. 

In [5] a so-called rejection-free version for general worm-type algorithms was 
proposed. If we consider a PS step in Al, then the total probability that one of 
the four conceivable moves takes place is given by 

A {u- k) = ^ + P2,+ + Pi,- + P2,-) . (41) 

Note that 1 — A > is the probability of keeping u unchanged. It is easy to show 
j5j that if we replace PS by a step PSnr 

• ii u = V holds, the same steps as in PS are taken, 

• for u ^ V one of the four moves discussed before is chosen and executed with 
the normalized probabilities 

then PSnr has the Boltzmann weight in 

Zn, = J2 t^-''<''"("''^'"V"' - v) 5[d*n - {5„>^ + (1 - 5^,,) A (u; k)} (43) 

n,u,v 

as a fixed point. To construct a complete simulation algorithm Alnr we need to 
replace beside PS — )■ PSnr also the CF step by one, where the Boltzmann ratio 



corresponding to (43 ) is formed for the acceptance step. Finally also a version A2nr 
is programmed, where A {u; k) is suitably generalized to depend on the direction 
of the defect lines. Since in any case the relative weights of vacuum graphs are the 



same in (27) and (43) all expectation values (23) before. 

We see that the term 'rejection free' refers to the Polyakov shifts in non-vacuum 
configurations only. Acceptance rates in CF are however at efficient levels in all 
our simulations. If now the remaining acceptance rate of PSnr out of vacuum gets 
too small, we may exploit our freedom of choosing p to improve this. For example 
the choice 

p{u - v) = 6ii^ij + c{l - 6^^^) (44) 

with c < 1 reduces the relative weight p^^ of vacuum configurations and thus 
enhances the probability to leave them. A few brief experiments allowed us to 
arrange for reasonable acceptance rates for these moves, too. We find this simple 
choice for p to be a good one at the critical point, where we will implement A2nr, 
while in the confined phase in connection with Al other choices will be more 
favorable. 
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4 Numerical Simulations 



4.1 Simulations close to the critical point 

We first report on a number of simulations at 

^, = 0.76141346(7) ^ tc = 0.64190876 (4) (45) 

wliicli corresponds to tlie dual of tlie estimate given for the spin model ii^ [22] . 
On each lattice of sizes L„ = L = 6, 8, 12, 16, 24, 32, 48 we have executed 8 x 10^ 



iterations of A2nr with p from (44) and c = (8/L) after a few tests. All observed 
acceptance rates of CF are around 50% and the remaining acceptance of PS out 
of vacuum configurations drops from 34% at L = 6 to 23% at L = 48 with these 
choices. 

One of the simplest observables is the average fraction of plaquettes partici- 
pating in the random surface of vacuum configurations 

^ = W- (46) 



which is analogous (and equivalent, see (25)) to the average plaquette. We measure 



values like for example 9 = 0.33751 (2) at L = 6 and 6 = 0.33493(4) at L = 
48. In Fig. [l] we visualize a typical graph occurring in the L = 6 simulation, 
larger systems look even more cluttered. Histograms counting the frequencies of 
separations u — v are accumulated after each individual PS step and turn out to 
be rather flat, as expected in the critical theory, where the Polyakov correlation 
decays only slowly. 

During each iteration of A2nr we separately accumulate contributions to ob- 
servables from the subset of vacuum configurations encountered during this itera- 
tion, like for example for B. In the end we perform an autocorrelation analysis of 
this time series to estimate errors as described in [21]. The such defined integrated 
autocorrelation time of is shown in Fig. |2} 

It is in units of A2nr iterations which cost proportional to L^. We see that while 
the absolute values are not too large on our lattices the rate of growth exhibits 
critical slowing down Tint oc with an effective dynamical exponent for our range 
of lattice sizes close to two. Note that we cannot make statements on the truly 
asymptotic dynamical exponent on the basis of these data. 



Of greater interest than G are the topological observables (18). Taking into 



account the symmetries between all planes we here measure the set of observables 

))o, n = 0,l,2,3. (47) 



^In the meantime the precision of the estimate for the critical temperature has been further 
improved |23| . 



15 



u-1 " 

Figure 1: A typical configuration on a 6^ lattice at t = tc- The Polyakov line 
defects are shown as fat red lines. 



Due to (24) each Rn can be simply related to ratios of partition functions with 
twisted and untwisted boundary conditions and is hence expected to possess finite 
continuum or scaling limits. In particular, at the exact critical point there are 
definite finite values i?* for each in the limit L oo aX t = tc- In fact, this 
property may be used to determine tc- Note that the values Rn are not indepen- 
dent, of course, but have to sum to unity. Our raw results on these quantities are 
compiled in Table [3] in the Appendix. We here consider the expected finite size 
scaling behavior close to the critical point in the form 

Rn {tc, L) = Rl + anL-^ + . . . , (48) 

where u is the exponent of the leading scaling violations. In our simulations at 
t = tc we find the results shown in Fig. [sj In the variable L~'^ a linear approach 
to Rn is thus expected. Recent determinations of the exponent u can be found in 
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data 
z=2.02 



10- 



12 16 24 32 

L 



48 



Figure 2: Double logarithmic plot of integrated autocorrelation times of G in units 
of A2nr iterations which are comparable to 'sweeps'. 



[22] and [25] - we use u = 0.82. 

We show simple fits for our lattices with L ^ 12. Details on raw data are 
given in the appendix. We here just list the extrapolated values Rq = 0.2306(10), 
Rl = 0.3754(6), R*2 = 0.3043(6), R^ = 0.0897(3). 



4.2 Simulations in the confined phase 

In the confined phase /3 < /3c we simulate with algorithm Al, such that the dis- 
tribution of u — V has a direct relation to the temporal Polyakov loop correlator 



(29). We do not generate surfaces wrapping around the 12 plane corresponding 
to dynamically summing over twisted and untwisted boundary conditions for this 
plane. The Polyakov correlator is however periodic in space in this situation. 

For our test we took the coupling f3 = 0.73107 with Li = L2 = L = 64 and a 
series of inverse temperatures Lq = 6, 8, ... , 26, 28. Such values (but smaller Lq) 
have been adopted in [T7| and a string tension in lattice units of a = 0.0440 (3) 
has been estimated based on earlier literature. In our simulations in the confined 



phase we want to use a p (x) in (27) that essentially captures the decay of the two 
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Figure 3: Plots of i?„ versus L ^ for simulation ai t = of [22] . 



point function (29). Inspired by the lowest term in (35) we prepar^ 



p (f ) oc ^ 



M 



Xi - kiLf + (x2 - k2Lf' 



P(0) 



(49) 



The mass M is estimated from a by using the first three terms on the right hand 



side of (34). The sum over integer which makes p periodic, is actually truncated 



after a few terms as the remainder would be insignificant compared to roundoff 
due to the exponential decay of the Bessel function. We emphasize that any 
imperfection in p is related to the efficiency of the simulations but does not entail 
any systematic error. 

For each value of Lq we have performed 1.2 x 10^ iterations of Al. In Fig. |4] 
we show the histogram of the frequencies of sampled separations of the two defect 



^ The sum is singular for xi — X2 — 0. We replace it by an extrapolation along the diagonal 
solving /(0, 0)//(l, 1) = /(I, l)//(2, 2). 



18 







Figure 4: Log of the number of occurrences of separations u — v of the Polyakov 
defect hues for the run 24 x 64^. 



hues in our run for Lq = 24. With the exception of very short torus distances 
we find a rather fiat behavior (note the fine scale of this logarithmic plot) which 
leads to very uniform bin heights, which, up to correlation effects, would naively 
imply constant relative errors. Via ([29]) this allows to measure the correlator which 
for the same run is plotted in Fig. [5] To determine errors here, we had to store 
the time series for these separations rather than just the histogram. We see that 
the exponential decay can be followed over the whole lattice with relative errors 
staying small. The plot jointly exhibits 'on axis' and 'diagonal' separations which 
form a very smooth curve and thus exhibit small violations of rotation invariance. 
The spatial correlation length is of order one and therefore the periodicity is just 
barely visible close to r = L/2 and r = L/ \/2. 

The closed string mass gap Eq is estimated from (36) and our rather pre- 
cise results are compiled in Table [T] A typical case of the determination of the 
mass gap is shown in Fig. |6} Effective masses m(y + 1/2) are determined by 
solving C (y) /C {y + 1) = cosh (m {y — L/2)) / cosh {m{y + 1 — L/2)). The hori- 
zontal lines show the error band from a fit deeply in the plateau region which leads 
to the entries for £"0 in the table. Observed autocorrelation times T^^t for these 
derived quantities range from 0.508(4) (Lq = 6) to 0.598(6) (Lq = 28). The 
remarkable feature here is, that the errors in the effective masses do not grow with 
separation. This is due to our judicious adaptation of the simulated ensemble by 
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Figure 5: Polyakov line correlation function (29) for the 24 x 64^ lattice versus 
Euclidean separation r = xf + x^. Circles refer to on-axis separations, crosses 
to diagonal ones. Since errors are invisible in the upper plot, we separately show 
the slowly growing relative errors in the lower plot. 



choosing p. The situation is analogous to the one in the Ising spin model ^ where 
the reasons for this success are discussed in more detail. 

What is the implication of our fluctuating boundary conditions in the 12 plane 
for the correlation (36)? The falloff in the direction xi is associated with the trans- 
fer matrix in this direction which acts on wave-'functionals' ip [c"o (^^o; ^2) , (a^o; 2:2)] • 
The Polyakov operator (28) acts on them by multiplication. The integration over 
links (Ti in the Euclidean theory implies the inclusion of a projector on gauge in- 
variant states in the thermal trace given by the path 'integral'. This refers to gauge 
transformations that are periodic in both directions. One may however in addition 
consider transformations [TH] that are topologically nontrivial, here antiperiodic in 
X2, and physical states may be even or odd under them. By summing over twisted 
and untwisted boundary conditions in the 12 plane we include a projector also 
with respect to this quantum number. Under the assumption that the ground 
state is even, such a projection is no disadvantage. 
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Lo 


Eo 


1^0 1 


Lo 


Eo 


1^0 1 


6 


0.160286(15) 


0.165454(28) 


18 


0.762774(77) 


0.104107(89) 


8 


0.279241(21) 


0.138969(31) 


20 


0.853968(98) 


0.10092(12) 


10 


0.384096(28) 


0.126905(37) 


22 


0.94450(12) 


0.09807(13) 


12 


0.482528(36) 


0.118760(45) 


24 


1.03480(14) 


0.09571(15) 


14 


0.577660(47) 


0.112776(55) 


26 


1.12416(16) 


0.09329(18) 


16 


0.670782(60) 


0.108026(60) 


28 


1.21429(19) 


0.09154(21) 



Table 1: Closed string mass gaps Eo and matrix element vo, see (31), on Lo x 64^ 
lattices at /3 = 0.73107. Each data point derives from 1.2 x 10'' iterations of Al. 



To analyze our data we perform a fit of the form 



Ei 



1 

u 



(50) 



To have visible errors at all we immediately plot the difference between the left 
and right hand side of such a fit against Lq ^ in Fig. [jj The fit here was derived 
from the subset Lo ^ 12 of the data. The point with Lq = 8 is much higher 
while Lq = 6 is far off the panel. At these time extents the low temperature 
expansion has clearly and rather abruptly broken down. In fact, for our /3 the 
phase transition is close to Lq ~ 4 according to [T7|. 



We note that (50) is exactly equivalent to the Nambu Goto form (32) if we 
identify 



Nambu Goto : 



TT 



Co 



a 



Cl 



-G 



_9_c| 

Co 



1. 



(51) 



In Table 2 we list a number of fits (50) with free co, ci together with the resulting 



ratio r for whose error the correlation between co and ci is taken into account. 



Lo,ram 


xV dgf 


r 




Cl 


10 


8.0/8 


1.0109(13) 


0.0440330(25) 


0.046363(35) 


12 


8.0/7 


1.0114(23) 


0.0440334(32) 


0.046374(62) 


14 


8.0/6 


1.0111(40) 


0.0440332(42) 


0.04637(11) 



Table 2: Fits (50) where Lo < Lo^ram are omitted. 



We see a small but significant deviation of r from one by only about 1%. The 
most likely explanation for this are cut off effects in our opinion, given that the 
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Figure 6: Effective mass in the run with Lq = 24. 



string tension a in lattice units, which may be talcen as indicative for their size, is 



about 0.044. It remains to be discussed to which terms in the expansion (33) our 



fits are actually sensitive. The errors of our Eq are 7%, 24% , 68% of the last term 
in (33) for Lq = 10,12,14. The analogous numbers for the next term are 50%, 
236%, 897%. This next term reads Stt^/ (10368(j^Lg) and this is only taken as a 
model of conceivable next order terms. We thus have to conclude that in spite of 
our rather high precision we are just still sensitive to the last (known) universal 
term of order Lq ^ and cannot confirm or disprove the coefficient of Lg here. 

In future simulations it will be desirable to simulate smaller lattice spacings to 
check the size of cutoff effects and possibly to enlarge the statistics to the point of 



seriously probing another term in (34). 



5 Conclusions 

We have attempted to generalize the worm algorithm for the Ising model to gauge 
theories with Ising link variables. The labeling of all-order (in tanh /3) strong cou- 
pling graphs as configurations of constrained plaquette fields was straight forward. 
The essential method of to achieve efficient updates of the graphs was to allow 
for a pair of point defects. The possible defects in the gauge case form a very 
large set of generalized loop networks. We could not identify among them a suit- 
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Figure 7: Difference between data of Table [T] and the second fit of the form ( 50 ) 
listed in Table [2l 



able subset that strongly reduces critical slowing down while staying 'close to' the 
vacuum or simple Wilson loop defects. 

The defects in the spin model allowed at the same time for a very precise 
estimation of the fundamental correlation at large distance [2]. We successfully 
have generalized this aspect to the gauge model and could compute the Polyakov 
line correlator with similar precision. At the same time, the algorithmic efficiency 
is rather high on large lattices in three dimensions in spite of critical slowing down. 

Other Abelian gauge theories have been considered in the surface representa- 
tion, see for example [7j, [9], ^U\. The techniques for correlators in this paper can 
clearly be generalized to these cases, i.e. to the groups Z(A^) or U(l). Theories 
with SU(A^) variables, both gauge and the principal chiral spin model^are open 
problems in this context to which we hope to come back in the future. 
Acknowledgments: UW is indebted to many people for helpful discussions: Tim 
Garoni, Martin Hasenbusch, Stefano Lottini, Mike Peardon, Stefan Sint, Erhard 
Seller, Rainer Sommer and Peter Weisz. Financial support of the DFG via SFB 
transregio 9 is acknowledged. 

® Note that while there is mention of worm simulations of SU(3) spin models in the literature 
|26j . their nontrivial symmetry is Z(3). 
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Appendix 



In this appendix we list the raw data on which Fig. |3]is based. Each hne represents 
8 X 10® iterations of algorithm A2nr as described in subsection 3.4 except for L = 48 



where the number has been doubled. The error of these quantities seems to grow 
roughly proportionally to L^'^ with a fixed number of iterations whose costs scale 
hke 



L 


i?0 


Ri 


R2 


Ri 


6 


0.23194(15) 


0.38290(10) 


0.29937(11) 


0.08579(6) 


8 


0.23112(19) 


0.38053(13) 


0.30103(13) 


0.08731(7) 


12 


0.23047(29) 


0.37920(18) 


0.30196(19) 


0.08838(10) 


16 


0.23075(40) 


0.37828(25) 


0.30227(26) 


0.08870(14) 


24 


0.22985(68) 


0.37813(41) 


0.30317(44) 


0.08884(23) 


32 


0.22983(110) 


0.37699(62) 


0.30338(69) 


0.08980(35) 


40 


0.23204(163) 


0.37553(90) 


0.30343(100) 


0.08901(50) 


48 


0.23182(165) 


0.37653(89) 


0.30321(100) 


0.08844(48) 



Table 3: Results for the observables (47) from simulations on L lattices with 



periodic boundary conditions in all directions at the estimated critical (3 given in 
(|45|. 
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